Probability distribution of the maximum of a smooth temporal signal 
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We present an approximate calculation for the distribution of the maximum of a smooth stationary 
temporal signal X{t). As an application, we compute the persistence exponent associated to the 
probability that the process remains below a non-zero level M. When X{t) is a Gaussian process, 
our results are expressed explicitly in terms of the two-time correlation function, f{t) = {X{0)X{t)). 



The problem of evaluating the distribution of the max- 
imum of a time-correlated random variable X(t)h.as 
elicited a large body of work bvmathematicians 111 1 3^ l3 , 
and physicists, both theorists B EEJl S 11 [llllirill 
and experimentalists [isl lisl lid Il7j |. In the physics 
literature, this is related to the persistence problem, the 
probability that a temporal signal X (and hence its max- 
imum) remains below a given level M up to time t. Per- 
sistence properties have been measured in as different 
systems as breath figures 0|, liquid crystals 0|, laser- 
polarized Xe gas 1 151. fl uctuating steps on a Si surface 
[ifi], or soap bubbles [lj|. 

The mathematical literature has mainly focused on 
evaluating P<(t) = Prob(X(t') < M, t' e [0,i]) for 
Gaussian processes and for large | M | , a regime where ef- 
ficient bounds or equivalent have been obtained jj, 13 . 
Recently 3], and for Gaussian processes only, a numeri- 
cal method to obtain valuable bounds has been extended 
to all values of M, although the required numerical effort 
can become quite considerable for large t. 

Physicists have also concentrated their attention to 
Gaussian processes 0, H ^ : which are often a good 
or exact description of actual physical processes. For in- 
stance, the total magnetization in a spin system |ll| . 
or the height profile of certain fluctuating interfaces 
0, 0, 0| are true temporal Gaussian processes. Two 
general methods have been developed, focusing on the 
case M = 0, which applies to many physical situations. 
The first one 0, Sill is a perturbation of the considered 
process around the Markovian Gaussian process, which 
has been extended for small values of M |^. Within 
this method, only the large time asymptotics of P<{t) 
is known, leading to the definition of the persistence ex- 
ponent (see below). The alternative method, using the 
independent interval approximation gives very ac- 
curate results for smooth processes, but is restricted to 
M = 0. 

In addition, this problem has obvious applications in 
many other applied and experimental sciences, where one 
has to deal with data analysis of complex statistical sig- 
nals. For instance, statistical bounds of noisy signals 
are extremely useful for image processing (for instance in 
medical imaging or astrophysics [Til), in order to obtain 



cleaner images by correcting spurious bright or dark pix- 
els 0,1^. In general, it is important to be able to evaluate 
the maximum of a correlated temporal or spatial signal 
originating from experimental noise. The same question 
can arise when the signal lives in a more abstract space. 
For instance, in the context of genetic cartography, sta- 
tistical methods to evaluate the maximum of a complex 
signal has been exploited to identify putative quantita- 
tive trait loci 0|. Finally, this same problem arises in 
econophysics or finance, where the probability for a gen- 
erally strongly correlated financial signal to remain below 
or above a certain level is of great concern. 

One considers a general {i.e. not necessarily Gaus- 
sian) stationary process X{t) of distribution g{X), and 
zero mean. The process is assumed to be "smooth", so 
that its velocity X'{t) is continuous, ensuring that the 
number of times N{t), where X = M, remains finite for 
any bounded time interval [0,t]. For a given level M, one 
defines r as the average temporal interval between two 
crossings of the level X = M. One also introduces P-{t) 
(respectively P-|_(t)) as the distribution of time intervals 
during which X{t) < M (resp. X{t) > M). The average 
of P± (t) are denoted by r± . Finally, one defines P< (t) 
(resp. P>{t)) as the probability, starting at X{0) < M 
(resp. ^(0) > M), that the process X remains below 
the level M (resp. above the level M) up to time t. The 
distribution of the maximum (resp. minimum) of the 
process X{t) in the interval [0,t] is clearly the derivative 
of P<(i) (resp. P>{t)) with respect to M. The difficulty 
of obtaining analytic forms for the above quantities lies 
in the fact that powerful methods like the Fokker-Planck 
approach are useless for non-Markovian processes. In this 
letter, we obtain closed expressions for P<(i), P>{t), and 
P± (t) from a minimal knowledge of the statistical prop- 
erties of the process X. This is achieved by analyzing 
the trajectories of X and using the sole assumption that 
the lengths of the intervals between successive crossings 
of the level M are uncorrelated. 

Let us assume that one knows the two following quan- 
tities A(t) and N^{t) from experiment, numerical simu- 
lations, or even analytically: A{t) is the autocorrelation 
function of 6[M — X(t)] (9 is Heaviside's function). 



A{t) = {6[M - X{t)] 0[M - X{Q)]), 



(1) 



2 



and 7V< (t) is the average number of crossings at level M 
up to time t, averaged over the starting position X{0) < 
M. For large time, one has 

iV<(t)~7V(t) = -, (2) 



with 



N{t) 



\X'{t')\S{X{t')~ M)dt' 



(3) 



If the process is smooth, X'{t') is not correlated 
with X{t'), and using stationarity, we find N{t) = 
tx{\X'{t)\){S{X{t) - M)), which leads to 



5(M)(|X'W|). 



(4) 



In addition, iV>(t), the average number of crossings at 
level M up to time t, starting from X{0) > AI, satisfies 
the sum rule 

G(M)iV<(t) + (1 - G{M))Ny{t) = -, (5) 

T 

where G{M) = g{x) dx. Note that t± are simply 
related to t = : 



2tG(M), 



2r(l - G{M)). 



(6) 



In the following, we obtain closed forms for P<{t), 
P>(i), and P±{t) from the knowledge of A{t) and iV<(t), 
for any level M. When X is a Gaussian process of corre- 
lator f{t), we shah see later that A(t), N^{t), Ny{t), r, 
and T± can be explicitly written in terms of /. Hence, the 
minimal knowledge of the two-time correlation function 
of a Gaussian process will grant access to the yet un- 
known quantities P<{t), P>{t), and P±{t). However, the 
present approach has a wider range of applications and 
does not rely on the Gaussian property of the process. 

Our central approximation consists in assuming that 
the interval length between crossings are uncorrelated 
[lot . The probability P<{N,t), starting from X{0) < M, 
that there are exactly N crossings in the interval [0,t], 
can then be written, for odd A^ = 2n— l(n>l), 

P<(2n-^l,i) =TrVo^*i<3-(^i)x 
/,* dt2 P+{t2 - ti) dt, P_ (i3 - t2) • • • 

jt2n.-2 "^^^n-l P- (^2n-l — t2n-2)Q +{i — i2n-l), (7) 

where Q±{t) = J^°° P±{t') dt' is the probability that a 
± interval is larger than t. For even N = 2n {n > 1), one 
obtains a similar expression 

P<{2n,t) ^ tZ^ dtiQ^{ti)x 

Jl dt2 P+{t2 - ii) !l dt3 P_ {t-s -t^)--- 

•i't2n-2 "^^Sn-l P-(i2n-l — ^2n-2)x 
Jt2„_i '^^^n P+{t2n — t2n-l)Q-{t — t27i)- (8) 



For any function of time F{t) introduced in this 
letter, one defines its Laplace transform F{s) = 
J^°° F{t) e""* dt. The convolution products in Eqs. Tfl^ 
take a much simpler form in the Laplace variable s 



P<(2n-l,s) = 
P<(2n,s) = 

_ 1-P±is) 



tZ'QIP+[P+P.] 



n-1 



(9) 
(10) 



One can now express the con- 



where Q±{s) 

servation of probability, P<(i) -I- J2~n^i P<i^^'^) 
which leads to 



In fact, Eq. Hll|l is an exact relation, which reads 

/ + CO 
{f -t)p_{t')dt', 



1, 



(11) 



(12) 



in the time variable. Indeed, if X{t) has not crossed the 
level M up to time i, it belongs to a — interval of duration 
t' > t, starting at an initial position uniformly distributed 
between and t' — t. Note that P>{s) and P>(iV, s) are 
given by similar expressions as Eqs. (|9I1UI11I) by exchang- 
ing the indices — and 

Now, P±is) can be calculated by expressing the known 
quantities A{s) and N^{s) as a function of P±(s): 



N<{-s) = 
A{s) = G(M) 



(1+P+)(1-P-) 
r_ s2(l - P+P_) 

1 1-^^+ 



I + P4 



.iV<(s) 



(13) 
(14) 



Using P4(0) — —T± and Eq. ©, one obtains the follow- 
ing estimates, valid for small s, 



N<is) 



A{s) 



(15) 



The first expression in Eq. H15|) is equivalent to Eq. lI2Jl, 
whereas the second relation expresses that for large 
A{t) ~ G2(M). For large s, 7V<(s) - [2G[M)ts^]~^ , 
which corresponds to the small time behavior 

t _ , , t 



N<{t) 



2G(M)t' 



N>{t) 



2(1 - G{M))n 



(16) 



For G{M) ^ \ {i.e. M ^ 0, when g{X) is symmetric), 
Eq. (|16|l differs from the large time asymptotics given by 
Eq. lI2Jl. However, the sum rule of Eq. Q is preserved by 
the small time estimates of Eq. (|16|l . Finally, writing 



F{s) 



G(M)^sA{s) 
G{M)sN<{s) ' 



(17) 



3 



and using Eqs. (|13I14() . the interval distributions read 



P+(s) 



l + FisY 

2-r-S^jV<(s)(l + F(s)) 
2 -r_s27V<(s) (1-/^(5))' 



(18) 
(19) 



Inserting these expressions of P± in Eq. (|ll|l , one obtains 
our final result for P< (and P>), from the sole knowledge 
of A{t) and N^{t) (or Ny{t)). 

The persistence exponent 9 is defined as the asymp- 
totic decay rate of P<(<) ~ e~**. The term "expo- 
nent" arises from the fact that in ma ny p hysical systems 
lllillSIIlIIlIIlIIlIIlIIlIIiEII, the process X 
of interest is stationary in the variable t = InT, where T 
is the actual physical time. Thus the persistence decays 
as a power law P<{T) ~ T~^, as a function of the real 
time T (see below for practical examples). Within our 
approach, —9 is the first pole of P<(s) (or equivalently 
of P_(s)) on the negative real axis 10]. Using Eq. H19[l. 
one finds that 9 satisfies the implicit equation 



9G{M)[1 + 9N<{' 



'-A{-9) 



(20) 



When M is large, 9 goes to zero, and Eq. H20|l leads to 

9 = tZ^ = {2tG{M))-\ (21) 

For a Gaussian process, the same expression was obtained 
from a heuristic argument in 3]. In this limit of large M , 
the interval distributions are found to become Poissonian. 
We conjecture that the present approach becomes exact 
for large M, the — intervals being so large that the -|- 
intervals are indeed uncorrelated. 

Let us move on to the case where X{t) is a stationary 
Gaussian process. The properties of X{t) (and hence 
A{t) and N^{t)) are completely determined by the sole 
knowledge of its two-time correlator f{t) = {X{0)X{t)). 
For a general process, this connection is only approxi- 
mate and can only be made by means of the IIA [ioj i. 
For convenience, we set {X'^{t)) = /(O) = 1. The process 
is smooth if / is twice differentiable. We also assume that 
for large time t, the correlator fit) decays fast enough so 
that f{t)dt is finite [l[l3. For t > 0, the position- 
velocity correlator is {X{Q)X' [t)) — f'{t), which vanishes 
for t — 0, since f{t) is an even function, twice differen- 
tiable at t = 0. The velocity-velocity correlation func- 
tion is {X'{0)X'(t)) = -f"{t). The mean time interval 
between crossings t is computed using Eq. Q 



vEZ!Me-^. 



(22) 



For a Gaussian process, A{t) has been derived in |2l| 



A{t)^ 



^ M - xf{t) 



-l{x)G 



dx. 



(23) 



For large time, so that f{t) is small, one finds 

Oifit)). (24) 



A{t)^G^{M) + l^e-''" 
2tt 



Finally, N^{t) can be calculated after introduc- 
ing the correlation matrix of the Gaussian vector 
{X{t),X{0),X'(t)), which reads 



1 fit) 
C{t) = \ fit) 1 fit) 
fit) -/"(O) 



(25) 



One finds 



iV<(t)=G-i(M) / {\X'it')\)^dt\ (26) 



where (|X'(t)|)< is the average of the velocity modulus, 
knowing that Xit) = M, and averaged over X(0) < M: 



i\X'it)\). 



M 



dxr 



— 00 — 00 



+00 



dv ■ 



me 2 



(27r)3/2^/d^' 



(27) 



where U = (M, xo,v). In practice, this integral has to be 
computed numerically, but can be reduced to a cumber- 
some one-dimensional integral over xq, involving g and 
G. For large time or small fit), one obtains 



{\X\t)\). 



1 



Mg(M) 
G(M)t 



fit) + Oifit)). (28) 



Note that for numerical purposes, the Laplace transform 
of it) can be efficiently written as 



N<is) 



1 

s .10 



+ 00 



i\x'it)\). 



1 



dt 



1 



(29) 



The various analytical asymptotic forms obtained above 
can be useful to complement the partial knowledge of 
Ait) and -/V<(i) from a partial experimental or numerical 
sampling of the process Xit) ■2(j\. 

As an application, Table reports the theoretical and 
numerical values (the latter obtained by direct simulation 
of the temporal process) of the persistence exponent 9 
for M — 0, 1, 2, 3 and for the two Gaussian processes 
Xi and X2 associated to the correlators 



hit) 

/2(i) 



1 



cosh 
1 



3c- 



(30) 
(31) 



Xi and X2 are two examples of non-Markovian Gaus- 
sian processes arising from physical systems. Indeed, 
up to a multiplicative term, Xiit) can be shown to be 
equal to /9(x, T) (at some arbitrary position x), where 
the density field p(x, T) evolves accor ding to the two- 
dimensional diffusion or heat equation [lOl |. ^ = V^/o, 



4 



0th Qsim Qsim 

0.1862 0.188(1) 0.2647 1/4 

1 5.914 ■ 10'^ 5.91(1) ■ 10"^ 8.625 ■ 10'^ 8.31(4) • 10"^ 

2 1.084-10"^ 1.09(1) -10"^ 1.665-10"^ 1.61(2) • 10~^ 

3 8.769 ■ 10"" 8.77(2) ■ 10"* 1.420 • 10"^ 1.42(1) • 10"^ 

TABLE L Exponents 9 from theory (S"*) and simulations 
(S^™), for different values of M, calculated for the processes 
X\ and X2 introduced in the text. For M = 0, the three 
first results were first reported in uA. while the exact value 
6»2(M = 0) = i was obtained in [gj"^ 

starting from an arbitrary (although not too strongly 
correlated) initial condition. Here, the actual time T 
is again related to our stationary time t by the relation 
t = InT. Note that the out of equilibrium dynamics of 
a two-dimensional Ising model after a quench at a tem- 
perature To < Tc (where Tc is the ferromagnetic criti- 
cal temperature) can be approximately mapped to this 
problem [loll2ll |. by assimilating the spin S'(x, T) = ±1 
to 5(x, T) = sign[/9(x,T)]. We find that the theoret- 
ical values are within the numerical error bars, except 
maybe for M = 0. Overall, the accuracy is better than 
1%. In addition, the asymptotic result of Eq. (|21|l al- 
ready leads to fair estimates for M — 2 and M = 3 
(6*1 (M = 2) = 1.102-10"2 and 6*1 (M = 3) = 8.852-10"'*). 
Our theoretical and numerical results are also consis- 
tent with the numerical bounds computed in 3], for 
M = 1 and M = 2 (0.0586 < ei{M = 1) < 'o.0684 
and 0.0106 < 0i(M = 2) < 0.0119). 

As for X2, it is associated with the random acceleration 
process j^^l, = rjiT), where rjlT) is a J-correlated 
white noise, in the variable T — e*. Contrary to Xi, the 
process X2 is not infinitely differentiable (/2"(i) is not 
defined at t — 0), although it is just smooth enough for 
the present approach to be applicable. Hence, it is not 
surprising that the theoretical results are not as good as 
for the process Xi. However, the theoretical estimates 
are clearly becoming more accurate as M increases, and 
presumably exact for large M. 

In conclusion, the present work develops a power- 
ful approximation leading to explicit expressions for the 
Laplace transform of the probability to remain above or 
below a certain level M (and hence the distribution of the 
minimum or maximum of the process). This approach 
also gives the distribution of time intervals during which 
the process remains above or below M, and leads to the 
determination of the persistence exponent. 

I am very grateful to Satya Majumdar and Partha Mi- 
tra for fruitful discussions. 
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